2 ノンパラ対応なし群間差
今回から群間の差の検定とそのグラフです。
一元配置のデータで、個体間要因とか級間要因とか betwee-factor とか言われる場合です。
まずはノンパラから。
2群の比較と3群以上の比較は統計のかけ方が異なるのですが、、、
グラフを描く上では2群もそれより多い場合も何も変わらないので、同じページで解説しています。
2群の比較の検定は、サンプルデータから2群だけ抽出する方法の解説とともにやります(結局、そのやり方をしないとできないので...)。
ノンパラは、平均値の比較ではななくて、数値データを順位に変換して、順位が高いデータがどの群に多いか?ということを比較していると思うので、平均値のグラフを使わない方が妥当だと思っています。なので、分布を可視化する表現方法の箱ひげ図(boxplot)を僕は使うことにしています。
箱ひげの、箱がその群の順位の1/4区切りを示していて、箱の内部に収まるデータがその群の半数で、真ん中の線が中央値です。上下の棒(これがヒゲ)の先端が最大値と最小値になることが多いですが、外れ値がある場合はヒゲの外にもdotsがあります。
https://scrapbox.io/files/638798a69072ae002237dd59.png
グラフの描画はそんなに難しくないのですが、実は、群間差の場合のノンパラ検定は、不等分散の場合に使えない(第1種の過誤がメチャ起る)ので、紹介はするものの、あまり使わないで頂きたい感じ。
code:r
# ---------- preparation in advanca 事前準備--------------------------------------------------
# Install packages if not installed 必要なパッケージが無けれな自動でインストール-----
if(!require("ggplot2")) install.packages("ggplot2")
if(!require("exactRankTests")) install.packages("exactRankTests") #needed fo Exact Rank Tests if(!require("lawstat")) install.packages("lawstat") #needed for Brunner-Munzel test # clear R's brain 一度参照したリストもしくは変数を次の作業で消す-----
rm(list = ls())
# Load libraries 必要なパッケージの読み込み-----
library(ggplot2)
library(NSM3) #needed for Steel-Dwass test library(exactRankTests) #needed fo Exact Rank Tests library(lawstat) #needed for Brunner-Munzel test # Data input 今後以下のパスからデータを自動で読み込みする場合-----
library(readr)
CanRsample <- read.csv(file.choose()) #Run and choose CanRsample2021.csv # 【2】box plots between 箱ひげ図 個体間----------------------------
kruskal.test(CanRsample$aveTemp_year ~ CanRsample$place) #(データシート$従属変数 ~ データシート$独立変数) pairwise.wilcox.test(CanRsample$aveTemp_year, CanRsample$place, paired=F, correct=F, exact=F, p.adj="holm")
#post-hoc multiple comparison, correct=F で連続性補正をしないようにしている、Rのデフォ設定でn50以上になると正規近似されるが、 pSDCFlig(CanRsample$aveTemp_year, as.factor(CanRsample$place), method = "Asymptotic")
dSendai <-filter(CanRsample, place == "Sendai") #(元データ, 抽出したい群名がある列名 == "群名") dTokyo <-filter(CanRsample, place == "Tokyo")
dKagoshima <-filter(CanRsample, place == "Kagoshima")
wilcox.test(dSendai$aveTemp_year, dKagoshima$aveTemp_year, paired=F, exact=F, correct=F)
wilcox.exact(dSendai$aveTemp_year, dKagoshima$aveTemp_year, paired=F)
ggplot(CanRsample, aes(x = place, y = aveTemp_year, group = place)) +
scale_x_discrete(limit=c('Sendai', 'Tokyo', 'Kagoshima')) +
geom_jitter(shape = 21, size = 4, aes(fill = place), width = 0.2, alpha = 0.7) +
scale_fill_manual(values = c (Sendai = "deepskyblue3", Tokyo = "gray1", Kagoshima = "orange" )) +
ggtitle("【2】box plots between") +
theme_classic(base_size = 22, base_family = "HiraKakuProN-W3") # white backgrond グラフの背景を白に
【統計の手法の確認】
パラメトリックな検定をする際には、
one-way ANOVA → Tukey-Kramer HSD 検定
という流れが用いられます。
ノンパラで3群以上の差を検定する際によく使われるのは、
Kruskal-Wallis test → Steel-Dwass test
という手順です。
Kruskal-Wallis test は、原理的には、マン・ホイットニーのU検定と同様だそうです。ちなみに、U検定とウィルコクソンの順位和検定(Wilcoxon rank sum test)は、同じ結果をもたらすものであると知られています。で、はじめにこのサイトについて の最後の方で書いた通り(と上でも行った通り)、不等分散に非常に弱い。 ただ、このページでは(かつての僕がせっかく調べたので)一通りやっていきます。
Steel-Dwassをしてくれる pSDCFlig という関数が入っている NSM3 というパッケージをこのページではインストールしています。僕が初めてNSM3を使った時は、さらに、partitions と gmp というパッケージがないと動かないとRに言われたのでインストールしたのですが、今回、これらが無くても動きました。一応、誰かのためにコメントとして残してあります。
さて、統計を進める手順は上記の通りなので、上から実行してもらえればOKです。
後検定の多重比較としては、Steel-Dwass検定以外にも、ウィルコクソンの順位和検定を総当たり(ペアワイズ)にやった上でHolm法でP値を調整をする方法も載せています。ウィルコクソンの用法(オプションの指定方法)も、スクリプト中のコメントを読んでもらえれば分かるかと思います。
【2群比較をする場合】
2群の比較の検定をする前に、dplyrパッケージのfilterを使って、群(Sendai、Tokyo、Kagoshima)ごとのデータシートのオブジェクトを作っています。なんでこんなことしてるのかと言うと、こういう指定の仕方をしないと、検定をしてくれない関数だからです。初めて使った時は、ちょっと信じられませんでした。2群だけが含まれるデータシートのオブジェクトだったとしても、ダメです。
検定ごと(関数ごと)に群・条件(独立変数)の指定の仕方が全然違うので、最初はとまどったいうか、調べるのに時間がかかったものです。
上記 script 中では、SendaiとKagoshimaを比較する場合の書き方だけ示しています。
Brunner-Munzel 検定についてだけ、最後に。
これは、U検定と違って、等分散性を仮定しなくて良い検定です。むっちゃ便利ですね。使える時は使うと良いと思います。ただ、たまに計算が収束しません。それと、3群以上の場合に使える似たような検定も、たぶんまだ知られてないと思います。そんなにメジャーじゃないので、生物系の雑誌で使った際に、Reviewerが知らない可能性があるので、ちょっと怖い。心理系では、だいぶ知られているのではないでしょうか。
【グラフ描画】
すでに 1_1 と 1_2 を読んでくださった方なら、簡単だと思うので、それら前のページには書いてない特徴だけ説明します。
今回は、boxplot() という関数で描きます(その名の通りですね)。
横軸の並び
散布図じゃなくて群間の差なので、横軸にはカテゴリカル変数(群の名前)が並びます。
その並びを自分の好みの順番に指定するための
scale_x_discrete()
というものを入れてあります。不要なら消してください。
dotの表示:散らすか一直線にするか
geom_jitter():dotを散らしたいときに使います(反復測定じゃないならこれが良いと思います、特に、nがでかいほどに)
geom_point():dotを散らさず一直線にするときに使います(反復測定で同一個体とかに線を引くならこっちが良い)
このどちらか1つを使います。
dotの色指定
グラフの場合、色は scale_fill_manual() で決めます。
scale_color_manual() ではない!という点に注意。